Untitled

1 Preamble

1.1 Library

pacman::p_load(tidyverse, janitor, sf, leaflet, leafem, stars, viridis)

mycity_loc <- read_csv("cities_in_malaysia.csv")

Rows: 22 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city
dbl (2): lat, lng

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mycity_pop <- read_csv("cities_population.csv") 

Rows: 22 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city
num (1): population

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mycity_comb <- left_join(mycity_loc, mycity_pop) %>% 
  mutate(shape = case_when(population < 500000 ~ 16, 
                           population >= 500000 ~ 15),
         colour = case_when(population > 200000 ~ "red", 
                            population >= 50000 & population <= 200000 ~ "orange",
                            population < 50000 ~ "yellow"))

Joining with `by = join_by(city)`

mycity_sf <- st_as_sf(mycity_comb, coords = c("lng", "lat"), crs = 4326)

mysf_nat <- st_read("gadm41_MYS_shp", layer = "gadm41_MYS_0")

Reading layer `gadm41_MYS_0' from data source 
  `C:\Users\puter\Documents\RStudio\24-06_KursusSpatial\gadm41_MYS_shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99.64072 ymin: 0.85372 xmax: 119.2697 ymax: 7.380556
Geodetic CRS:  WGS 84

mysf_state <- st_read("gadm41_MYS_shp", layer = "gadm41_MYS_1")

Reading layer `gadm41_MYS_1' from data source 
  `C:\Users\puter\Documents\RStudio\24-06_KursusSpatial\gadm41_MYS_shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99.64072 ymin: 0.85372 xmax: 119.2697 ymax: 7.380556
Geodetic CRS:  WGS 84

mysf_dist <- st_read("gadm41_MYS_shp", layer = "gadm41_MYS_2")

Reading layer `gadm41_MYS_2' from data source 
  `C:\Users\puter\Documents\RStudio\24-06_KursusSpatial\gadm41_MYS_shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 144 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99.64072 ymin: 0.85372 xmax: 119.2697 ymax: 7.380556
Geodetic CRS:  WGS 84

2 Visualization

2.1 MY City Location

mycity_loc %>% 
  ggplot(aes(x = lng, y = lat)) +
  geom_point() + 
  labs(x = "Longitude", y = "Latitude", title = "City Location in Malaysia") + 
  theme_bw()

mycity_loc %>% 
  ggplot(aes(x = lng, y = lat)) +
  geom_point(colour = "red", 
             size = 5, # base plot = pch
             shape = 16) + # base plot cex
  labs(x = "Longitude", y = "Latitude", title = "City Location in Malaysia") + 
  theme_bw()

2.2 MY City Population

mycity_comb %>% 
  ggplot(aes(x = lng, y = lat, 
             colour = colour, shape = as.character(shape))) + 
  geom_point(size = 5) + 
  scale_color_manual(values = c("red" = "red", 
                                "orange" = "orange", 
                                "yellow" = "yellow")) + 
  scale_shape_manual(values = c("15" = 15, "16" = 16)) +
  labs(x = "Longitude", y = "Latitude", title = "City Population in Malaysia") + 
  theme_bw()

2.3 Adjust Bounding Box

if use ggplot, there is no need to use bounding box function to adjust the bounding box.

ggplot(data = mycity_sf) +
  geom_sf(aes(colour = colour, shape = as.factor(shape)), size = 5) +
  coord_sf(xlim = c(90, 120), ylim = c(-5, 20)) +
  scale_color_manual(values = c("red" = "red", 
                                "orange" = "orange", 
                                "yellow" = "yellow")) +
  scale_shape_manual(values = c("16" = 16, "15" = 15)) +
  labs(x = "Longitude", y = "Latitude", title = "City Population in Malaysia") +
  theme_bw()

2.4 Add MY Shapefile

ggplot(data = mycity_sf) +
  geom_sf(data = mysf_state, fill = "grey", colour = "black") +
  geom_sf(aes(colour = colour, shape = as.factor(shape)), size = 5) +
  scale_color_manual(values = c("red" = "red", 
                                "orange" = "orange", 
                                "yellow" = "yellow")) +
  scale_shape_manual(values = c("16" = 16, "15" = 15)) +
  labs(x = "Longitude", y = "Latitude", title = "City Population in Malaysia") +
  theme_bw()

2.5 Leaflet Package

leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data = mysf_state,  
              color = "green", 
              weight = 1, 
              fillOpacity = 0.2, 
              popup = ~NAME_1) %>% 
  addCircleMarkers(data = mycity_comb, 
                   lng = ~lng, lat = ~lat, 
                   radius = ~sqrt(shape), 
                   color = ~colour, 
                   fillOpacity = 0.8, 
                   popup = ~city)

2.6 stars Package

2.6.1 World

world_elevmap <- read_stars("wc2.1_2.5m_elev.tif") %>% 
  as.data.frame(., as_points = F, na.rm = T)

ggplot(world_elevmap) +
  geom_raster(aes(x = x, y = y, fill = wc2.1_2.5m_elev.tif)) +
  scale_fill_viridis_c() +  # Optional: Enhance color scale
  labs(title = "Elevation World", fill = "Elevation") +  # Add titles and labels
  theme_minimal()  # Use a minimal theme

2.6.2 Malaysia

mal_elevmapst <- read_stars("wc2.1_2.5m_elev.tif") %>% 
  st_crop(., mysf_nat)

Warning in st_crop.stars(., mysf_nat): st_crop: bounding boxes of x and y do
not overlap

mal_elevmapdf <- mal_elevmapst %>% 
  as.data.frame(., as_points = F, na.rm = T)

mal_elevmapdf %>% 
  ggplot() + 
  geom_raster(aes(x = x, y = y, fill = wc2.1_2.5m_elev.tif)) +
  scale_fill_viridis_c() +  # Optional: Enhance color scale
  labs(title = "Elevation World", fill = "Elevation") +  # Add titles and labels
  theme_minimal()  # Use a minimal theme

2.7 Leaflet and Stars

elev_values <- as.vector(mal_elevmapst[[1]])  # Extract the raster values
color_palette <- colorNumeric(palette = viridis(256, option = "C"), 
                              domain = range(elev_values, na.rm = TRUE))


leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addStarsImage(mal_elevmapst, 
                band = "wc2.1_2.5m_elev.tif", 
                project = F, 
                opacity = 0.8,
                colors = color_palette) %>% 
  addPolygons(data = mysf_state,  
              color = "green", 
              weight = 1, 
              fillOpacity = 0.2, 
              popup = ~NAME_1) %>% 
  addCircleMarkers(data = mycity_comb, 
                   lng = ~lng, lat = ~lat, 
                   radius = ~sqrt(shape), 
                   color = ~colour, 
                   fillOpacity = 0.8, 
                   popup = ~city)